import java.util.Arrays;
import java.util.Random;

public class zuoye1 {
    private int nComponents;
    private double[] mean;
    private double[][] components;
    private double[] explainedVarianceRatio;

    public zuoye1(int nComponents) {
        this.nComponents = nComponents;
    }

    public void fit(double[][] X) {
        System.out.println("=== PCA手动实现：三维数据降维到二维 ===");

        // 1. 数据标准化（去中心化）
        this.mean = calculateMean(X);
        double[][] XNormalized = normalizeData(X, mean);

        System.out.println("原始数据:");
        printMatrix(X);
        System.out.println("均值: " + Arrays.toString(mean));
        System.out.println("标准化后的数据:");
        printMatrix(XNormalized);

        // 2. 计算协方差矩阵
        double[][] covarianceMatrix = calculateCovariance(XNormalized);
        System.out.println("协方差矩阵:");
        printMatrix(covarianceMatrix);

        // 3. 计算特征值和特征向量
        EigenResult eigenResult = eigenDecomposition(covarianceMatrix);
        double[] eigenvalues = eigenResult.eigenvalues;
        double[][] eigenvectors = eigenResult.eigenvectors;

        System.out.println("特征值:");
        for (int i = 0; i < eigenvalues.length; i++) {
            System.out.printf("λ%d = %.6f\n", i + 1, eigenvalues[i]);
        }

        System.out.println("特征向量:");
        for (int i = 0; i < eigenvectors.length; i++) {
            System.out.printf("v%d = [", i + 1);
            for (int j = 0; j < eigenvectors[i].length; j++) {
                System.out.printf("%.6f", eigenvectors[i][j]);
                if (j < eigenvectors[i].length - 1) System.out.print(", ");
            }
            System.out.println("]");
        }

        // 4. 选择主成分
        this.components = selectComponents(eigenvalues, eigenvectors);
        System.out.println("选择的主成分:");
        printMatrix(components);

        // 5. 计算解释方差比例
        this.explainedVarianceRatio = calculateVarianceRatio(eigenvalues);

        System.out.println("解释方差比例:");
        for (int i = 0; i < explainedVarianceRatio.length; i++) {
            System.out.printf("主成分 %d: %.4f (%.2f%%)\n",
                    i + 1, explainedVarianceRatio[i], explainedVarianceRatio[i] * 100);
        }

        double totalVariance = 0;
        for (int i = 0; i < nComponents; i++) {
            totalVariance += explainedVarianceRatio[i];
        }
        System.out.printf("前%d个主成分保留的方差: %.4f (%.2f%%)\n",
                nComponents, totalVariance, totalVariance * 100);
    }

    public double[][] transform(double[][] X) {
        if (mean == null || components == null) {
            throw new IllegalStateException("必须先调用fit方法训练模型");
        }

        double[][] XNormalized = normalizeData(X, mean);
        return matrixMultiply(XNormalized, transpose(components));
    }

    public double[][] fitTransform(double[][] X) {
        fit(X);
        return transform(X);
    }

    private double[] calculateMean(double[][] X) {
        int nSamples = X.length;
        int nFeatures = X[0].length;
        double[] mean = new double[nFeatures];

        for (int j = 0; j < nFeatures; j++) {
            double sum = 0;
            for (int i = 0; i < nSamples; i++) {
                sum += X[i][j];
            }
            mean[j] = sum / nSamples;
        }
        return mean;
    }

    private double[][] normalizeData(double[][] X, double[] mean) {
        int nSamples = X.length;
        int nFeatures = X[0].length;
        double[][] normalized = new double[nSamples][nFeatures];

        for (int i = 0; i < nSamples; i++) {
            for (int j = 0; j < nFeatures; j++) {
                normalized[i][j] = X[i][j] - mean[j];
            }
        }
        return normalized;
    }

    private double[][] calculateCovariance(double[][] X) {
        int nSamples = X.length;
        int nFeatures = X[0].length;
        double[][] covariance = new double[nFeatures][nFeatures];

        for (int i = 0; i < nFeatures; i++) {
            for (int j = 0; j < nFeatures; j++) {
                double cov = 0;
                for (int k = 0; k < nSamples; k++) {
                    cov += X[k][i] * X[k][j];
                }
                covariance[i][j] = cov / (nSamples - 1);
            }
        }
        return covariance;
    }

    private EigenResult eigenDecomposition(double[][] matrix) {
        int n = matrix.length;
        double[] eigenvalues = new double[n];
        double[][] eigenvectors = new double[n][n];

        // 保存原始矩阵用于特征值计算
        double[][] originalMatrix = copyMatrix(matrix);
        // 复制矩阵用于收缩
        double[][] currentMatrix = copyMatrix(matrix);

        Random random = new Random(42); // 固定随机种子以便重现

        for (int comp = 0; comp < n; comp++) {
            System.out.printf("  正在计算第 %d 个特征向量...\n", comp + 1);

            // 初始化随机向量
            double[] eigenvector = new double[n];
            for (int i = 0; i < n; i++) {
                eigenvector[i] = random.nextDouble();
            }
            eigenvector = normalize(eigenvector);

            for (int iteration = 0; iteration < 1000; iteration++) {
                // 矩阵乘以向量
                double[] newVector = matrixVectorMultiply(currentMatrix, eigenvector);

                // 归一化
                newVector = normalize(newVector);

                // 检查收敛（只比较方向，不比较大小）
                if (distance(newVector, eigenvector) < 1e-10) {
                    break;
                }

                eigenvector = newVector;
            }

            // 使用原始矩阵计算特征值（瑞利商）
            double eigenvalue = rayleighQuotient(originalMatrix, eigenvector);

            // 存储结果
            eigenvalues[comp] = eigenvalue;
            for (int i = 0; i < n; i++) {
                eigenvectors[comp][i] = eigenvector[i];
            }

            System.out.printf("    找到特征值: %.6f\n", eigenvalue);

            // 矩阵收缩（Hotelling收缩法）- 使用当前矩阵
            if (comp < n - 1) {
                double[][] outerProduct = outerProduct(eigenvector, eigenvector);
                currentMatrix = matrixSubtract(currentMatrix,
                        scalarMatrixMultiply(eigenvalue, outerProduct));
            }
        }

        // 按特征值大小排序
        Integer[] sortedIndices = argsort(eigenvalues, false); // 降序
        double[] eigenvaluesSorted = new double[n];
        double[][] eigenvectorsSorted = new double[n][n];

        for (int i = 0; i < n; i++) {
            eigenvaluesSorted[i] = eigenvalues[sortedIndices[i]];
            eigenvectorsSorted[i] = eigenvectors[sortedIndices[i]];
        }

        return new EigenResult(eigenvaluesSorted, eigenvectorsSorted);
    }

    private double[][] selectComponents(double[] eigenvalues, double[][] eigenvectors) {
        // 按特征值降序排列
        Integer[] sortedIndices = argsort(eigenvalues, false);
        double[][] selectedComponents = new double[nComponents][eigenvectors[0].length];

        for (int i = 0; i < nComponents; i++) {
            selectedComponents[i] = eigenvectors[sortedIndices[i]];
        }
        return selectedComponents;
    }

    private double[] calculateVarianceRatio(double[] eigenvalues) {
        double totalVariance = 0;
        for (double eigenvalue : eigenvalues) {
            totalVariance += eigenvalue;
        }

        double[] varianceRatio = new double[eigenvalues.length];
        for (int i = 0; i < eigenvalues.length; i++) {
            varianceRatio[i] = eigenvalues[i] / totalVariance;
        }
        return varianceRatio;
    }

    // ========== 数学工具方法 ==========

    private double[] normalize(double[] vector) {
        double norm = norm(vector);
        if (norm < 1e-12) return vector; // 避免除零
        double[] normalized = new double[vector.length];
        for (int i = 0; i < vector.length; i++) {
            normalized[i] = vector[i] / norm;
        }
        return normalized;
    }

    private double norm(double[] vector) {
        double sum = 0;
        for (double v : vector) {
            sum += v * v;
        }
        return Math.sqrt(sum);
    }

    private double distance(double[] a, double[] b) {
        double sum = 0;
        for (int i = 0; i < a.length; i++) {
            double diff = a[i] - b[i];
            sum += diff * diff;
        }
        return Math.sqrt(sum);
    }

    private double[] matrixVectorMultiply(double[][] matrix, double[] vector) {
        int n = matrix.length;
        double[] result = new double[n];

        for (int i = 0; i < n; i++) {
            for (int j = 0; j < n; j++) {
                result[i] += matrix[i][j] * vector[j];
            }
        }
        return result;
    }

    private double rayleighQuotient(double[][] matrix, double[] vector) {
        double[] mv = matrixVectorMultiply(matrix, vector);
        double numerator = dotProduct(vector, mv);
        double denominator = dotProduct(vector, vector);
        return numerator / denominator;
    }

    private double dotProduct(double[] a, double[] b) {
        double result = 0;
        for (int i = 0; i < a.length; i++) {
            result += a[i] * b[i];
        }
        return result;
    }

    private double[][] outerProduct(double[] a, double[] b) {
        int n = a.length;
        double[][] result = new double[n][n];

        for (int i = 0; i < n; i++) {
            for (int j = 0; j < n; j++) {
                result[i][j] = a[i] * b[j];
            }
        }
        return result;
    }

    private double[][] matrixSubtract(double[][] a, double[][] b) {
        int rows = a.length;
        int cols = a[0].length;
        double[][] result = new double[rows][cols];

        for (int i = 0; i < rows; i++) {
            for (int j = 0; j < cols; j++) {
                result[i][j] = a[i][j] - b[i][j];
            }
        }
        return result;
    }

    private double[][] scalarMatrixMultiply(double scalar, double[][] matrix) {
        int rows = matrix.length;
        int cols = matrix[0].length;
        double[][] result = new double[rows][cols];

        for (int i = 0; i < rows; i++) {
            for (int j = 0; j < cols; j++) {
                result[i][j] = scalar * matrix[i][j];
            }
        }
        return result;
    }

    private double[][] matrixMultiply(double[][] a, double[][] b) {
        int aRows = a.length;
        int aCols = a[0].length;
        int bCols = b[0].length;
        double[][] result = new double[aRows][bCols];

        for (int i = 0; i < aRows; i++) {
            for (int j = 0; j < bCols; j++) {
                for (int k = 0; k < aCols; k++) {
                    result[i][j] += a[i][k] * b[k][j];
                }
            }
        }
        return result;
    }

    private double[][] transpose(double[][] matrix) {
        int rows = matrix.length;
        int cols = matrix[0].length;
        double[][] result = new double[cols][rows];

        for (int i = 0; i < rows; i++) {
            for (int j = 0; j < cols; j++) {
                result[j][i] = matrix[i][j];
            }
        }
        return result;
    }

    private Integer[] argsort(double[] array, boolean ascending) {
        Integer[] indices = new Integer[array.length];
        for (int i = 0; i < array.length; i++) {
            indices[i] = i;
        }

        Arrays.sort(indices, (i, j) -> {
            if (ascending) {
                return Double.compare(array[i], array[j]);
            } else {
                return Double.compare(array[j], array[i]);
            }
        });

        return indices;
    }

    private double[][] copyMatrix(double[][] matrix) {
        double[][] copy = new double[matrix.length][matrix[0].length];
        for (int i = 0; i < matrix.length; i++) {
            System.arraycopy(matrix[i], 0, copy[i], 0, matrix[i].length);
        }
        return copy;
    }

    private void printMatrix(double[][] matrix) {
        for (double[] row : matrix) {
            System.out.println(Arrays.toString(row));
        }
    }

    // 特征值分解结果容器
    private static class EigenResult {
        double[] eigenvalues;
        double[][] eigenvectors;

        EigenResult(double[] eigenvalues, double[][] eigenvectors) {
            this.eigenvalues = eigenvalues;
            this.eigenvectors = eigenvectors;
        }
    }

    // ========== 主函数 ==========
    public static void main(String[] args) {
        // 固定的三维数据
        double[][] data3D = {
                {2.5, 2.4, 1.8},
                {0.5, 0.7, 1.2},
                {2.2, 2.9, 2.1},
                {1.9, 2.2, 1.9},
                {3.1, 3.0, 2.8},
                {2.3, 2.7, 2.3},
                {2.0, 1.6, 1.5},
                {1.0, 1.1, 0.9},
                {1.5, 1.6, 1.4},
                {1.1, 0.9, 1.0}
        };

        System.out.println("=".repeat(60));
        System.out.println("PCA手动实现：三维数据降维到二维");
        System.out.println("=".repeat(60));

        // 创建PCA模型
        zuoye1 pca = new zuoye1(2);

        // 执行PCA
        double[][] data2D = pca.fitTransform(data3D);

        System.out.println("\n" + "=".repeat(60));
        System.out.println("最终结果");
        System.out.println("=".repeat(60));
        System.out.println("降维后的二维数据:");
        pca.printMatrix(data2D);

        System.out.printf("\n原始数据形状: %d x %d\n", data3D.length, data3D[0].length);
        System.out.printf("降维后形状: %d x %d\n", data2D.length, data2D[0].length);

        // 计算信息保留率
        double totalVariance = 0;
        for (int i = 0; i < 2; i++) {
            totalVariance += pca.explainedVarianceRatio[i];
        }
        System.out.printf("信息保留率: %.4f (%.2f%%)\n", totalVariance, totalVariance * 100);

        // 验证关键计算
        System.out.println("\n=== 关键计算验证 ===");
        verifyKeyCalculations(data3D);
    }

    private static void verifyKeyCalculations(double[][] data3D) {
        // 验证均值计算
        double[] mean = new double[3];
        for (int j = 0; j < 3; j++) {
            double sum = 0;
            for (double[] sample : data3D) {
                sum += sample[j];
            }
            mean[j] = sum / data3D.length;
        }
        System.out.println("验证均值: " + Arrays.toString(mean));

        // 验证第一个样本的中心化
        double[] firstSample = data3D[0];
        double[] normalized = {
                firstSample[0] - mean[0],
                firstSample[1] - mean[1],
                firstSample[2] - mean[2]
        };
        System.out.println("验证第一个样本中心化: " + Arrays.toString(normalized));
    }
}